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SUMMARY 


The present work used an existing Navier-Stokes code (PARC2D) to 
compute the nozzle flow field. Grids were generated by the interactive 
grid generator codes TBGG and GENIE. All computations were made on the 
NASA/M5FC CRAY X-MP computer Comparisons were made between the 
computations and MSEC in-house wall pressure measurements for C0 2 

flow through a conical nozzle having an area ratio of 40. Satisfactory 
agreements exist between the computations and measurements for 
different stagnation pressures of 29.4, 14.7 and 74 psia, at stagnation 
temperature of 1060 °R. However, agreements did not match precisely 
near the nozzle exit. Several reasons for the lack of agreement are 
possible. The computational code assumes a constant gas gamma, whereas 
the gamma i.e. the specific heat ratio foi C0 2 varied from 1.22 in the 

plenum chamber to 1.38 at the nozzle exit. The computations also assumes 
adiabatic and no-slip walls. Both assumptions may not be correct. 
Finally, it is possible that condensation occurs during the nozzle expansion 
at the low stagnation pressure. The next phase of the work will 
incorporate variable gamrna and slip wall boundary conditions in the 
computational code and develope a more accurate computer code. 
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I. IN TRODUCTION 


During the last decade, the requirements for more accurate ana economical 

computation methods for predicting fluid flows in practical engineering 
problems has become increasingly demanding. Since full scale tests or 
model tests are often expensive and in many cases inconclusive, it is 
indispensable to rely on the analysis of computational fluid dynamics. It is 
also important to have adequate numerical methods and computer codes 
available to help the interpretation of experimental data and to aid vehicle 
design. The treatments of fluid flows have advanced tremendously in 
recent years due to the continuous improvements of advanced computer 
such as CDC, CYBER and the CRAY. 

In space applications the rocket engines for attitude controls and orbital 
maneuvers are often required to produce a thrust less than one pound force. 
The requirements of such low thrust engine dictate the use of low chamber 
pressures, low mass flow rates and small nozzles. It results that the 
nozzle flow is in the low Reynolds number range with viscous effects 
dissipated across the whole nozzle. An accurate knowleoge of th»ir thrust 
is required for designing the spacecraft control system. Recently the 
man-made contaminations to the spacecraft, from the nozzle exhaust has 
become an issue, and accurate solutions for the nozzle/plume are required. 
The flow in such nozzles possesses strong viscous/inviscio interactions at 
their exit due to the thick boundary layers. Traditional nozzle design 
techniques, such as the use of the method of characteristics to calculate 
the inviscid core and boundary layer theorey to compute the displacement 
thickness, fail to predict the strength of the viscous/inviscid interaction. 
Therefore, it is necessary to use a full Navier-Stokes code for this purpose. 
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Based on these premises, this project was proposed to develop an adequate 
numerical method of solvidng the full Navier -Stokes equations in three 
years ano to predict the two-dimensional mixed sjbsonic/supersonic flow 
characteristics of the low thrust converging/diverging nozzle, which the 
viscous effects are eminent across the entire nozzle. 

The present work used a full Navier -Stokes code (PARC2D) to compute the 
nozzle flow field Grids were generated using the interactive grid 
generator codes TBGG and GEN ,r All computations were made on the 
NASA/MSFC CRAY X-MP computer Comparisons were made between the 
computations and wall pressure measurements for cOj flow through a 

conical nozzle having an area ratio of 40 Satisfactory agreements exist 
between the computations and measurements for different stagnation 
pressures of 29.4, 14 7 and 7.4 psia, at stagnation temperature of 1060 °R. 
However, agreement did not match precisely near the nozzle exit. Several 
reasons for the lack, of agreement are possible. The computational code 
assumes a constant gas gamma, whereas the gamma i.e. the specific heat 
ratio for COj varied from 122 in the plenum chamber to 1.38 at the nozzle 

exit. The computations also assumes adiabatic and no-slip walls. Both of 
which may not be correct Finally, it is possible that condensation occurs 
during the nozzle expansion for low stagnation pressure The next phase of 
the work will incorporate variable gamma and slip wall boundary 
conditions in the computational code. Our first year's experience on 
running this existing code enhances the capability to develope a more 
accurate and efficient computer code for this type of flow problem 
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The first objective or this project is to compute the CO 2 nozzle flow fields 

by using existing computer code (PARC2D) and to compare the results with 
the experimental data. The second objective is to develop a computer code 
for solving the compressible full Navier-Stokes equations for low thrust 
type nozzle flow calculations The third objective is to include the 
rarefield gas dynamic in this code by adding the velocity and temperature 
slip conditions along the nozzle wall 

3 APPROACH 

3 1 Experimental Test Results 

A series of low thrust nozzle tests were performed by Or Lynn C. Chou and 
James A Carter [1] at NASA/MSFC in 1953. The experiment was on a CO 2 

nozzle flow field and its associated plume expanding into a highly rarefied 
space environment. The pressures along the nozzle wall and pitot 
pressures out of the nozzle along the center line will be used in this 
project as the reference data base. Figure 1 shows the nozzle 
configuration and the static pressure stations 

3 2 Remote Computer Setup 

A remote computer network was established under this project. An IBM 
personal computer system 2, model 50 with a modem (Fastalk) was choosen 
as remote terminal to connect to the EADS (Engineering And Design System) 
of super computer CRAY-X/MP AT NASA/MSFC The terminal emulator VTEK 
allows IBM PS/2 to emulate DEC VTI00 terminal as well as Tektromc 
4105/4014/4010 terminals. Figure 2 shows the arrangement through the 
communication network PSCN. 
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Figure 1 NOZZLE CONFIGURATION 
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Figure 2. REMOTE COMPUTER WORK STATION SETUP 
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3.3 Grid Generation 

The first step to calculating the flow in a nozzle is to generate an 
appropriate computational grid. Two grid generator codes have been used. 
The first grid generator is called TBGG and was developed by Smith and 
Wiese 12] at Langley Research Center. This grid generator was run on the 
VAX System in NASA/M5FC. It is interactive in the Tektronic 4010 mode 
and can be used to generate body fitted grids algebracally. The output file 
then was copieo to a magnetic tape and read onto the EADS System. The^ 
second grid generator is called GENIE and was developed by Soni {3] . GENIE 
code can provide computational grids for wide range of geometries related 
to internal flow problems. The process uses several techniques either 
separately or in combination to generate grids quickly and economically 
for arbitrary geometries. GENIE code was adopted later in this project due 
to it's more generic structure. 

Computations were performed on the nozzle configuration, the major 
geometric parameters for the nozzle and upstream plenum chamber are 
shown in Figure 3. This geometry was introduced into the grid generator as 
a discrete number of dimensional parts or as a set of equations. The 
exponential packing function was used to cluster the grids near the walls 
as shown in Figure 4. The blow up of grids near the nozzle throat region is 
shown in Figure 5. 

3.4 Governing Equations 

Two-dimensional Navier-Stokes equations in fully conservative form are 
used to describe the nozzle flow The transformed equations can be 
written as 

“r + + F„ + H ■ <R e + s„ + T) 
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Where 


u - u/j 

E - (C t U + E + F)/J 
F • ( n t u + n x E + Hy F)/J 
H - (F + H)/(yJ) 

R m <C X R + C y S)/J 
s • (n x r + n y s)/J 
T * (S + T)/(yJ) 


And the curvilinear coordinate system ^ andT are represented as 


t - t, i - c (t,x,y) , n - n (c.x.y) 
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The quantities r xx ,T X y.Tyy are components of the stress tensors, fts the 
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ratio of specific heat, k is the coefficient of thermal conductivity, Pr is 
the Prandtl Number and Re is the Reynolds Number. In order to complete the 
set of governing equations, the equation of state 
P-f R T 

is used, where ^ is the density, R is the gas constant and T is the absolute 
temperature. 

3.5 Navier-Stokes Solver 

- X. 

The Navier-Stokes code, PARC2D, is a modification of the ARC2D code that 
was developed by Pulliam and Steger [4] at NASA/ Ames Research Center. It 
uses a thin layer approximation to the Navier-Stokes equations with 
parabolized viscous stress terms. It is written in the strong conservative 
form in curvilinear coordinates The non-iterative implicit approximate 
factorization scheme of Beam and Warming [5],[6] is used with a fourth 
order artif ical dissipation. 

The PARC2D code is a modification of the ARC2D code by Cooper [7] at 
Sverdrup Technology, Inc./AEDC Group which removed the thin layer 
approximation. It is fully elliptic, requiring closed boundary conditions and 
an initial condition everywhere in the flow field. The code is modular and 
fully vectorized. It assumes that the gas is perfect gas with constant 
gamma and constant Prandtl number, and the Sutherland viscosity law for 
the temperature variation of viscosity. 
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Figure 3. GEOMETRY OF CONICAL NOZZLE AND PLENUM CHAMBER 
(ALL DIMENSIONS ARE IN INCHES) 
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FIGURE: 5 BLOW-UP OF GRID 90X81 GRID NEAR NOZZLE THROAT REGION 
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3.6 Initial and Boundary Conditions 

For the present calculations the stagnation conditions were taken as 
reference and the reference length was one inch, since the grid coordinates 
were calculated in inches. The Reynolds number is based upon a reference 
sound speed and length. As suggested by Collins 18} that the fastest 
convergence occurred by takinq u - v - o and density and energy equal to 
their stagnation values. These values were assigned initially throughout 
the flow field. 

Other options used are as follows: adiabatic wall, axisymmetric, viscous, 
laminar. Free boundaries (allowing inflow or outflow) were assumed at the 
upstream and downstream boundaries, and the other boundaries were axis 
of symmetry and no si ip/adiabatic wall. 

Boundary conditions must be imposed upon all of the boundaries. For the 
free boundaries this consists cf a specification of the pressure. The 
stagnation pressure was imposed upstream but the downstream boundary 
condition posed a problem. For flows with strong viscous/inviscid 
interaction the pressure cannot be expected to be constant across nozzle 
and greater pressure variations would exist across the exit of a conical 
nozzle. However, a constant downstream pressure must be specified. It is 
usually prudent to specify a pressure somewhat lower than the minimum 
expected downstream pressure. No boundary conditions are required on the 
axis nor on the wall. If constant temperature wail conditions are assumed 
then the temperature must be specified. Note that the entire wall does not 
have to be assumed to be at the same temperature. 

An efficient convergence procedure for those nozzle problems was also 

suggested by Collins. The steps are as follows: 
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1) Initially set q - 0, /> * 1 and E - 1/rtK- 1) everywhere in flow field. 

2) Set initial parameters as follows: DIS2 - 0.2, DIS4 - 0.35, PCQMAX - 
10.0, DTCAP - 5.0. Run until the axial velocity is positive everywhere 
in the plenum chamber. 

3) Slowly reduce DIS4 to 0.25 (all other parameters constant). 

4) Slowly reducde DIS2 to 0.00 

5) Check DT and set DTCAP to about one-half of the minimum DT for the 
last series of interations. Then run until L2 reaches an acceptalbe 
value (1 0" 8 to 10~ 9 ). 

In this discussion DIS4 and DIS2 are parameters related to the fourth order 
and second order dissipation, respectively, PCQMAX sets the maximum 
change in any variable during an iteration, DT is the time step and DTCAP is 
the maximum allowable time step. L2 is a convergence measure. 

4. RESULTS 

All the computations were made on the CRAY-X/MP super computer at 
NA5A/MSFC. A grid of 90x81 was used in the calculations which required 
approximately 800k computer storage. Convergence occurred after 7000 
iterations. Computations were made using the test conditions that 
corresponded to the MSEC in-house measurements on the nozzle described 
In Figure l The test gas was CO 2 , stagnation temperature was 1060*R and 
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the stagnation pressure had three values, 29.4, 14.7 and 7.4 psia. Since 
only wall pressure measurements were made, that was the only parameter 
that could be used for comparison with the computations. The 

computational code assumed a constant gamma whereas the gamma for C0 2 

varied from about 1.22 in the plenum chamber to 1.38 at the exit as shown 
in Figure 6. Therefore, exact comparison could never be expected to occur 
between the measurements and computations. 

The calculated results was illustrated by the PLOT3D, which is a 
commercial software package to plot the flow field distributions. PLOT3D 
is available at NASA/MSFC which can be run on either EADS or IRIS 
workstation. The flow parameters, such as pressure, temperature, velocity 
and Mach number can be plotted as contours on the nozzle crossection. In 
order to illustrate some particular parameters at a specific location (such 
as the Mach number along the center line), some specific plotting programs 
were developed by using graphical software package DISSPLA to fit these 
needs. 

Comparisons of computed and measured wall pressure at different chamber 
nressures are given in Figures 7, 8 and 9 respectively. The agreements at 
each stagnation pressure are quite well. It is difficult to distinguish the 
calculated and measured results from these figures. However, it shows 
some differences near the nozzle exit in the semi-log scale in Figure 10, 
which is drawn to magnify the differences. Several reasons are possible 
for the lack of agreement, in addition to the need for a variable gamma 
computational code. These include the possibility that the adiabatic wait 
boundary condition is not applicable and computations at a constant wall 
temperature should be performed to examine that possibility. In addition, 
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there is a strong possibility that condensation of the C0 2 was occurring 

during the expansion and agreement would not be expected until the 
stagnation temperature was raised to eliminate all possibisity for 
condensation. Finally, at the exit the Knudsen number is about 0.06, based 
on a mean free path at wall conditions and a crude estimate of the 
displacement thickness. Slip is expected to occur under such 
circumstances and the slip wall boundary condition formulated by Collins 
[91 will be implimented . 

Typical profiles of Mach number and static temperature across the nozzle 
at a location near the nozzle exit are shown in Figures 1 1 and 12. The large 
temperature gradients near the wall are caused by the assumed adiabatic 
wall condition. 

Figure i3 shows the Mach number along the center line. The Mach number 
increases rapidly af tei passed the throat region. The contour plots of Mach 
number, Pitot Pressure, Temperature and Velocity in Figures 14, 15, 16 and 

17 are presented to show the computational features of the entire field. 
These contour plots clearly show the phenomena of strong viscous/ inviscid 
interaction along the nozzle wall, especially near the nozzle exit. Figure 

18 shows the velocity-vector field on the upper half of the nozzle. It is 
seen that the boundary layer thickness is increasing along the nozzle flow 
stream. These results show that the present code predicts the low thrust 
nozzle quite well except near the nozzle exit. 
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Figure 6. THE SPECIFIC HEAT RATIO OF C0 2 GAS 
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Figure 7. WALL PRESSURE VS CENTER DISTANCE 
(CHAMBER TEMPERATURE- 1 060°R, 
CHAMBER PRESSURE *29 4 psia) 
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F igure 8. WALL PRESSURE VS CENTER DISTANCE 
(CHAMBER TEMPERATURE* ! 060°R, 
CHAMBER PRESSURE* 1 4 7 psia) 
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Figure 9. WALL PRESSURE VS CENTER DISTANCE 
(CHAMBER TEMPERATURE- 1 060°R, 
CHAMBER PRESSURES psia) 
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Figure 10 WALL PRESSURE VS CENTER DISTANCE 
(CHAMBER TEMPERATURE' 1 060°R) 
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Figure 1 1. MACH NUMBER VS RADIAL DISTANCE 
(CHAMBER TEMPERATURE* 1 060°R I 
CHAMBER PRESSURE- 29.4 psia) 
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Figure 12. TEMPERATURE VS RADIAL DISTANCE 
(CHAMBER TEMPERATURE * 1 060°R, 
CHAMBER PRESSURE- 29 A psia) 
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Figure 1 3. MACH NUMBER VS CENTER DISTANCE 
(CHAMBER TEMPERATURE - 1 060°R, 
CHAMBER PRESSURE- 29.4 psia) 
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Figure 14. MACH NUMBER CONTOUR LINES 
(CHAMBER TEMPERATURE- 1 060°R i 
CHAMBER PRESSURE- 29 4 psia) 
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Figure 15. PITOT PRESSURE CONTOUR LINES 
(CHAMBER TEMPERATURE- 1 060°R, 
CHAMBER PRESSURE- 29 4 psia) 
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Figure 16. TEMPERATURE CONTOUR LINES 
(CHAMBER TEMPERATURE - 1 060°R, 
CHAMBER PRESSURE* 29.4 psia) 


25 





-0.7S -0.50 -0.25 0.00 0.25 0.50 0.75 1.00 1.25 


1 


i 

\ 



Figure 17. HORIZONTAL VELOCITY CONTOUR LINES 
(CHAMBER TEMPERATURE* 1 060°R, 
CHAMBER PRESSURE* 29.4 psia) 
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Figure U. VELOCITY VECTOR 
(CHAMBER TEMPERATURE = 1 060°R, 
CHAMBER PRESSURE 3 29 4 psia) 
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5. CONCLUSIONS AND FUTURE WORKS 


The following conclusions have been made from the results of this study 

1) The grid generator code GENIE provided a good means generating useful 
computational grids for the Navier-Stokes computations of the nozzle 
flow 

2) The PARC2D code yields reasonable solutions for the flow field in 
supersonic nozzles at low Reynolds numbers where large 
viscous/inviscid interaction exists. Because the code is modular it can 
be easily modified 

A number of works will be done in the second year : 

1 ) The variable gamma to the PARC2D code will be implemented. 

2) The code will include the slip boundary conditions. 

3) In addition to the conical nozzle, the calculation will be performed on 

the sonic nozzle flow field and the plume, which is currently under 
testing in AEDC. The configuration is shown in Figure 19 and the 
computational grid has been developed as shown in Figure 20 This will 
be one of the tasks in the next year. 
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Figure 19 CONFIGURATION OF SONIC NOZZLE 



Figure 20. 90x61 GRID APPLIED TO THE SONIC NOZZLE 
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